Construction of analytical many body wave functions for correlated bosons in a 

harmonic trap 
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We develop an analytical many-body wave function to accurately describe the crossover of a one- 
dimensional bosonic system from weak to strong interactions in a harmonic trap. The explicit wave 
function, which is based on the exact two-body states, consists of symmetric multiple products of 
the corresponding parabolic cylinder functions, and respects the analytically known limits of zero 
and infinite repulsion for arbitrary number of particles. For intermediate interaction strengths we 
demonstrate, that the energies, as well as the reduced densities of first and second order, are in 
excellent agreement with large scale numerical calculations. 
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Introduction Ultracold dilute quantum gases repre- 
sent an amazingly rich platform for the realization of 
strongly interacting many-body systems [H. Extensive 
control of the trapping geometry via external fields as 
well as tuning of the interaction properties of ultracold 
atomic ensembles is nowadays routinely possible [2I, l3|. 
Theoretical models of quantum many-body systems that 
were considered as rough idealizations in the past, can 
now be prepared in a very pure way in order to study, 
for example, quantum phase transitions [4]. One very 
fundamental model of this type is explored in the present 
work: an ensemble of bosons at zero temperature, con- 
fined to a one-dimensional (ID) harmonic trap and in- 
teracting via a repulsive contact potential of arbitrary 
strength. Although conceptually simple, this system is 
in general not analytically solvable. However, the ex- 
perimental preparation of the trap is possible and the 
interaction strength can be tuned to any desired value 
by adjusting the strength of magnetic fields in the vicin- 
ity of a Feshbach resonance Q , or by employing confine- 
ment induced resonances a tool specific to quasi-lD 
waveguide- like systems. 

In view of the rich experimental possibilities, many of 
the ground-breaking theoretical works performed in the 
1960's, such as the celebrated Gross-Pitaevskii equation 
[gI], have attracted large attention. Focusing on ID sys- 
tems, and interaction strengths beyond the mean-field 
regime, two seminal works of this period are of special 
interest: (i) the Lieb-Liniger solution [7] via the Bethe 
Ansatz for contact interacting bosons in absence of ex- 
ternal potential with periodic boundary conditions (ex- 
tendable to hard- wall boundaries 0) and (ii) the Tonks- 
Girardeau gas [9] for hard-core bosons mapped to non- 
interacting fermions (Bose-Fermi map) through the ef- 
fective Pauli exclusion principle imposed by the infinite 
repulsion (fermionization). The Tonks- Girardeau gas [lO| 
and its counterpart for the attractive excited state, the 
Super- Tonks gas [11], were both experimentally realized 
within the last decade. The Bose-Fermi map works for 
arbitrary potential geometry, while the Lieb-Liniger so- 



lution applies to arbitrary interaction strength; in the 
homogeneous case, the two coincide for infinite repulsion. 

In this letter we construct an analytical many-body 
wave function which describes contact interacting bosons 
in a parabolic external potential for arbitrary interaction 
strength. It reproduces the two limiting cases of zero 
interaction and infinite repulsion for the harmonic trap. 
Our approach is based on products of functions inspired 
by the exact solution of the underlying two-body prob- 
lem. It shows an impressive agreement with results from 
extensive numerical studies and offers the possibility of 
extensions to, e.g., higher dimensions or other trap ge- 
ometries. 

Hamiltonian For the investigation of a ID trapping 
geometry one should take into account the experimen- 
tal conditions and their effect on the collision proper- 
ties of the atoms. Experimentally, the standard method 
to create quasi- ID tubes is using a very strong laser 
field for the transversal direction compared to the lat- 
eral one [l3, [ll|. This way the trap becomes highly 
anisotropic with the characteristic transversal length 

scale a_\_ = much smaller than the longitudi- 

[ijO_\_ (cj||) is the transversal (longitudi- 
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nal ay = 

nal) harmonic confinement frequency]. Then the trans- 
verse degrees of freedom are energetically frozen to their 
ground states and the effective ID interaction strength 

reads gio = (l - ^75^)"' where ao is the 

3D s-wave scattering length. 

The ID N-body Hamiltonian reads: 



H 



i<j 



where the contact interaction of particles located at 
Xi, i = 1,...,A^ is represented by the Dirac ^-function, 
and lengths and energies are scaled by ay and huj\\ re- 
spectively. The single remaining parameter is thus the 
rescaled interaction strength g = gijjy/AT/Wuj^^ which is 
altered either by tuning ao via magnetic Feshbach reso- 
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nances, or a± by modifying the transversal confinement. 

An important property of H that we employ here is the 
separability H = Hcm + Hr, where Hcm = ^ + 
describes a harmonic oscillation in the center-of-mass co- 

ordinates = ^ ~ ^ — ' while Hr 

i=i i=i 
describes the relative motion of the particles. The rel- 
ative motion is non-trivially coupled for any choice of 
relative coordinates for > 2, and will be exclusively 
addressed in the following 

Correlated pair wave function The correlated pair 
wave function (CPWF) developed here, is inspired by 
the idea that the pairwise contact interaction may be 
adequately addressed in the many-body system, if the 
discontinuity that it causes is imposed on each pair of 
atoms in the ensemble, in a similar way as for a single 
pair. The CPWF for the relative motion is formed as 
a pairwise product expansion of functions based on the 
two-body solutions, thereby respecting the two exactly 
solvable limits of zero and infinite interaction strength 
for an arbitrary number of particles. Specifically, the 
construction principle is composed of the following three 
postulates: 

(i) The CPWF for the relative motion of N particles is 
a product of parabolic cylinder functions (PCF) Df^ [13] 
of the distance rij = \xi — Xj \ of each pair: 
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Figure 1: (a) Wave function of the relative motion for two 
particles ^cp(r = xi — X2) and several strengths of the in- 
teraction (/i = 0, 0.5, 0.8, 1) and comparison of /i = 1 with 
the fermionic antisymmetric function. Contour plots of the 
density distribution |^(xi, X2, xs)!^ = 0.01 for three particles 
using the CPWF for the relative part for (b) /a = 0.2, (c) 
fi = 0.5, (d) fi = 0.8. 



^c^ = Cl[D^{l3r 



(1) 



i<j 



where P = ^^^""'"^ is the number of distinct pairs and the 
parameters /3 and /i, to be determined next, are identical 
for every pair since we deal with identical particles, while 
the absolute value enforces the bosonic permutation sym- 
metry. The normalization constant C will, without loss 
of generality, be omitted in the following, 
(ii) We make the assumption 
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which ensures that reproduces the exact analytically 
known solutions in the limits g = and g 00, for 
any number of particles. In these limits, /i equals and 
1 respectively (this follows from the boundary condition 
imposed next), and in both cases (as well as for every 

integer /i) D^{x) = e~~He^(x) where He^(x) are the 
modified Hermite polynomials. Therefore the total wave 

function ^ = ^i?^r with = e and = ^ 



reproduces the noninteracting ^ = exp 



N 2 ' 

— X7 



cp 



and 



N 



infinitely repulsive limit ^ = exp 



i=l / i<j 



[13|. The latter coincides with the (fermionic) determi- 
nant, written in form of pairs, and with implemented 
bosonic symmetry (Tonks- Girardeau limit). 

(iii) The boundary condition 2pD[.{^) = gDij{0), 

where D[j d{/3r%) D^{Prij), is imposed 

for each pair in the ensemble at r^j = 0, and determines 
fi for a certain value of the interaction strength g. With 

R. 
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the known expressions for the PCF I)^(x = 0) — — 



and 0) 
dental equation 



^ ^ ^ [12\ , the resulting transcen- 



2tr(V) 



(3) 



is solved for /i, selecting the solution in the interval ji G 
[0, 1] which corresponds to the ground state. The origin 
of this boundary condition will be discussed below. 

Discussion Postulate (i) addressing the construction 
of a many-body wave function via a product of functions 
of the relative distance appears also in other treatments 
in similar form (see Ref. "14)). Nevertheless, the particular 
choice of the PCF as a building block proves to be, as we 
will show, a very efficient approach for the present prob- 
lem, owing its inspiration to the two-body exact solution 
[l5| . Key features of our approach are postulates (ii) 
and (iii) determining the properties of the PCF, thereby 
ensuring that ^cp reproduces the two analytical limits 
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of zero and infinite coupling for any A^, and constitutes 
the exact solution at arbitrary interaction strength for 
N = 2. Furthermore, in the two-particle case the ana- 
lytical solution derived by Busch et. al [15] possesses a 
high degree of generality: it holds for bosons or fermions 
in one-, two- and three-dimensional harmonic traps and 
arbitrarily strong attractive or repulsive contact interac- 
tions, including ground and excited states of the relative 
motion. This readily implies that since our CPWF with 
this solution as a building block is sufficiently accurate 
for the ID repulsive bosonic case, a similar Ansatz could 
be envisaged for other setups. 

To clarify the implied boundary condition in postulate 
(iii), as well as why the CPWF is exact for two particles 
we investigate the action of the relative motion Hamilto- 
nian operator 



on 



d 

dxi 
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dXn 
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Introducing the notation 



ikl,mn,... 



Y[D^{f3rij^ki,mn,...) and xim = 1 - 2(^zm, we obtain: 



i<j 



Hr^cp = Y.^{nj){-2m^+gDii)4>'^ (4) 



P 



ikl,mn 



i<j 



k,ljm,n 



where sums over non-repeated terms with k ^ I ^ 

k,l,m,n 

m ^ n. The ffist sum in Eq. @ which contains the S- 
interaction for each pair vanishes at r^j = by imposing 
the boundary conditions 2PD[j{0) = gDij{0). This pro- 
vides the behaviour at the two-body contact point arising 
from the discontinuity of D[- for each pair. Therefore the 



manifolds A4 = {(xi, ...x^, Xj, xat) G J? 



ATI 



of contact (two particle collision) being of dimensionality 
N — 1 are taken into account correctly, while higher or- 
der contact (three particle collision etc.) represent lower 
dimensional manifolds. The second sum in Eq. dU is in 
the form of the Weber equation [l2[: D"[x) — ^D{x) = 
-(/i + l)D{x), for which the PCFs are solutions, and 
contributes to the energy by P(/i + ^). For N = 2 
(/3 = 1, P = 1) the last two terms in Eq. (j4]) vanish 
and therefore with the corresponding boundary con- 
dition is the exact solution with energy of the relative 
motion Cr = (/i+ ^). For N > 2 the last two sums in Eq. 
dH) do not vanish and provide additional contributions to 
the energy. 
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Figure 2: (a) Energies as a function of the interaction strength 
g obtained numerically and analytically via the CPWF. One- 
body densities pi(xi), for (b) three and (c) four particles for 
interaction strengths g = 0.2, 2.0, 10.0 comparing numerical 
and analytical results (only the xi > part is shown since it 
is symmetric with respect to the ordinate). In (b) a compar- 
ison with a variational calculation for intermediate ^ = 2.0 
is shown, (d) Two-body density p2{xi,X2) for three particles 
analytically calculated for interaction strength g — 20.0. 



Before comparing our analytical approach with cor- 
responding numerical results, we illustrate in Fig. [1] 
the spatial distribution of two and three particles. For 
two particles the wave function of the relative motion 
^cp(^ = xi — X2) acquires a cusp at r = for ^ > 
[or /i > from Eq. (j3j)] which goes to zero as ^ ^ 00 
(/i 1), retrieving the fermionic state but with bosonic 
symmetry. Similarly, for three particles the contour plots 
demonstrate, for increasing interaction strength, deple- 
tion of the probability density along the collision mani- 
folds {vij = 0); the conceptually important physical in- 
sight offered by our CPWF approach is that it captures 
the correlation properties in the vicinity of collision sur- 
faces, i.e., the tendency of the particles to repel and thus 
avoid each other. 

Accuracy of There are several numerical ap- 

proaches to the problem which are limited in terms of 
the particle number they can handle for a given accuracy: 
Exact diagonalisation [16|, Quantum Monte Carlo jlTj . 
and multi-configurational approaches in terms of Hartree 
products have been employed. We will resort here 
to the latter method which is reliable for a small num- 
ber of particles, assuring high numerical accuracy with a 
dense grid, small width of the Gaussian extrapolation of 
the ^-function and large basis set (we use 500 basis func- 
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tions and 30 single-particle orbitals). For details of the 
computation of stationary states by improved relaxation 
(propagation in the imaginary time) we refer the reader 
to Ref. 18 and references therein. 

The first observable we present in Fig. [2] (a) for three 
and four particles is the expectation value of the energy 
with increasing interaction strength g. We observe ex- 
cellent agreement between the energy calculated using 
CPWF and the numerical results, typically amounting 
to a relative accuracy of the order of 10 ~^ for interme- 
diate values of g. Remarkably, for larger g -a regime 
where convergence of numerical methods or approxima- 
tive models is challenging [19]- the analytical calculation 
yields lower ground state energies, and thus proves to be 
more accurate. 

Next we analyze reduced density opera- 
tors, namely one- and two- body densities, 
Pi{xi) J^^... J^^\'^\'^dx2...dxN and P2(^i,^2) 
J^^...J^^\^\^dx3...dxN. The CPWF is shown to 
capture very well not only qualitatively, but also quan- 
titatively, the properties of the crossover from weak 
to strong correlations. As a result of the increasing 
repulsion between the bosons, the one-body density [Fig. 
[2] (b),(c)] flattens and forms an A/'-peak structure (in- 
cluding xi < 0) close to fermionization |l3[ ll6Ml9| . Note 
that for intermediate interaction {g = 2.0) the deviation 
between CPWF and numerical results is somewhat 
larger; this is to be expected since is exact in, and 
therefore very accurate close to, the limits ^ = and 
g ^ oo {fi = and /i ^ 1), while between these limits 
(/i ^ 0.5) the error of the analytical approach should 
be maximal. Even further improvement is obtained 
by relaxing postulate (ii) and performing a variational 

optimization of f3 instead of keeping it fixed at 

[Eq. (|2])]. For optimal fSy^r ^ 0.84 the energy is slightly 
lower by 0.1% and the one-body density fits excellently 
to the numerical one [see Fig. [2] (b)]. Properties of the 
two-body density such its depletion along the diagonal 
xi = X2, due to the avoided simultaneous contact (a 
behavior attributed to the cusps of the CPWF in Fig. 
[T]), and peaks on the off-diagonal are shown in Fig. [2] 
(d) for three strongly interacting particles in agreement 
with the results in Refs. |l3|, |l8|. 

Conclusions The many-body wave function pre- 
sented in this letter for the problem of zero-range repul- 
sively interacting bosons with arbitrary coupling strength 
in a one-dimensional harmonic trap, has been inspired by 
exact analytical solution of the two-body case. We have 
shown that taking a product expansion of pair wave func- 
tions similar to the two-body solution, but modified such 
that the two analytically known limits of zero and infi- 
nite interaction strength are exactly reproduced for any 
number of particles, leads to an impressive agreement 
with extensive numerical calculations for any value of 
the couphng strength. Therefore the construction prin- 



ciple of our approach may be applicable to other setups, 
e.g. in higher dimensions, where the building block -the 
two-particle function- can be represented by an analytical 
expression similar to the exact solution of the relative mo- 
tion for the two-body problem, and some limiting cases 
of the relevant parameters (here e.g. zero and infinite 
interaction strength) are either exactly or approximately 
known. Improvements of our approach may be possi- 
ble (e.g. by taking into account lower dimensional spaces 
where three or more particles meet), and it is an excellent 
starting point for variational calculations (performed also 
in an exemplary case here) or for Quantum Monte-Carlo 
methods as an appropriate guiding function of Jastrow 
type, or for exploring the limits of mean- field treatments. 
Along with an accurate description of a highly correlated 
many-body problem, we believe that our approach offers 
valuable physical insight into the correlation properties 
at collision manifolds, which is conceptually useful in the 
theoretical treatment of many-body physics. 
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